Exercise 1 : Manipulate sf objects and associated data.frames
Exercise 2 : Static maps
Exercise 3 : Interactive maps

The main data used in this tutorial and in the lecture are about the geolocalisation of french restaurants in Paris. Data is extracted from an official register called SIRENE (Computer system for the business and establishment register) managed by the French National Institute of Statistics and Economic Studies (Insee) and geolocated by Etalab (French task force for Open Data). This register records the civil status of all companies and their establishments (including restaurants). SIRENE has the advantages of being rigorous and exhaustive on the French territory.

Exercise 1 : Manipulate sf objects and associated data.frames

1

Import the iris¹ map layer ‘iris_75.shp’ of Paris.

Clues

Use sf::st_read().

Solution

library(sf)
iris_75 <- st_read("data/iris_75.shp")

Reading layer `iris_75' from data source `/home/tim/Documents/prz/satRday/exercises/data/iris_75.shp' using driver `ESRI Shapefile'
Simple feature collection with 992 features and 4 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 643075.6 ymin: 6857477 xmax: 661086.2 ymax: 6867081
epsg (SRID):    NA
proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

2

Display the basemap of Paris with plot(iris_75). What do you notice ?

Solution

plot(iris_75)

Solution

We notice that R performs 4 graphs: one graph per variable in the sf object.

3

What is the functionality of the sf::st_geometry() function? What solution do you propose then?

Solution

sf::st_geometry() makes it possible to isolate the information contained in the ‘geometry’ column of the sf object. Using it, we put aside other variables (here CODE_IRIS, P14_POP, AREA and CODE_COM).

Solution

plot(st_geometry(iris_75))

4

Import the restaurant layer ‘sir.shp’ and display a map of Paris with its restaurants.

Solution

Use st_read() and sf::st_geometry().

Solution

sir_75 <- st_read("data/sir_75.shp")

Reading layer `sir_75' from data source `/home/tim/Documents/prz/satRday/exercises/data/sir_75.shp' using driver `ESRI Shapefile'
Simple feature collection with 23348 features and 7 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 643502.1 ymin: 6857645 xmax: 659766.4 ymax: 6867041
epsg (SRID):    NA
proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

plot(st_geometry(iris_75), bg = "cornsilk", col = "lightblue", 
     border = "white", lwd = .5)
plot(st_geometry(sir_75), col = "red", pch = 20, cex = .2, add=TRUE)
title("Restaurants in Paris")

5

Count the number of restaurant by iris.

Clues

Use sf::st_intersects() and sapply().

Solution

wwww

Solution

inter <- st_intersects(x = iris_75, y = sir_75)
iris_75$RESTAU <- sapply(inter, length)
head(iris_75)

Simple feature collection with 6 features and 5 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 650181.1 ymin: 6861761 xmax: 652179 ymax: 6863138
epsg (SRID):    NA
proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
  CODE_IRIS P14_POP      AREA CODE_COM                       geometry RESTAU
1 751010101     974  64006.07    75101 MULTIPOLYGON (((652096.5 68...     44
2 751010102     167  62671.06    75101 MULTIPOLYGON (((652052.7 68...      6
3 751010103     246  47653.41    75101 MULTIPOLYGON (((651869 6862...     16
4 751010104       3 222656.99    75101 MULTIPOLYGON (((651639.9 68...     10
5 751010105       0 243255.23    75101 MULTIPOLYGON (((650369.8 68...      0
6 751010199       0 234277.01    75101 MULTIPOLYGON (((652096.5 68...      0

6

Using the layer called ‘iris_75’, create a new aggregated map layer called ‘com_75’ which corresponds to the Paris Arrondissements. Also keep in this new layer the information on the population, area and number of restaurant in each municipality.

Information

The map layer called ‘iris_75’ contains the 5 digit codes of arrondissemnt in its variable CODE_COM.

Clues

Use the classic functions of dplyr package: select, group_by et summarize. These functions also work with sf objects.

Solution

library(dplyr)
com_75 <- iris_75 %>%
  group_by(CODE_COM) %>%
  summarize(P14_POP = sum(P14_POP), 
            AREA = sum(AREA), 
            RESTAU = sum(RESTAU))

Exercise 2 : Static maps

We would like here to design EPCI maps that combine the number of restaurants and the number of restaurants per 10,000 inhabitants. The EPCI correspond to the level of french intercommunalites.

1

Data preparation:

Load ‘fra.shp’ in R: the map layer of the borders of the French metropolitan territory.
Load the layer called ‘epci_31.rds’ (which contains the population and the number of restaurants in each EPCI) and create a variable called nb_rest_10000inhab which corresponds to the number of restaurants per 10,000 inhabitants in each territory.
Create a vector of quantiles breaks of the nb_rest_10000inhab variable.
Create the vector colors which corresponds to a the number of classes defined earlier.
Finally, add a variable called typo to ‘epci_31’ which indicates the class of the territory according to the discretization contained in bks for the nb_rest_10000inhab variable (this is needed for ggplot2 maps).

Information

For the creation of ‘bks’ et ‘cols’, use the getBreaks et carto.pal functions of the cartography package. For the creation of the typo variable, you can use the cut function and apply the parameters digit.lab = 2 and include.lowest = TRUE.

Solution

library(sf)
library(cartography)
library(dplyr)
# Import data
fra <- st_read("data/fra.shp", quiet = TRUE)
epci_31 <- readRDS("data/epci_31.rds")
# Create the variable
epci_31$nb_rest_10000inhab <- 10000 * epci_31$nb_of_rest / epci_31$P14_POP
# Define breaks
bks <- getBreaks(v = epci_31$nb_rest_10000inhab, method = "quantile", nclass = 4)
# Define color palette
cols <- carto.pal("orange.pal", length(bks)-1)
# Create a "typo"" variable
epci_31 <- epci_31 %>%
  mutate(typo = cut(nb_rest_10000inhab,breaks = bks, dig.lab = 2,
                    include.lowest = TRUE))

2

With the help of cartography package, make the following map which contains in a choropleth layer the variable nb_rest_10000inhab and in a proportional circle layer the variable nb_of_rest. The you can try to do the same map using the ggplot2 and tmaps packages.

With cartography:

Solution

# Define plot margins
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
# Find EPCI bounding box
bb <- st_bbox(epci_31)
# Plot France using EPCI boundingbox
plot(st_geometry(fra), col="ivory", border = "ivory3",
     xlim = bb[c(1, 3)], ylim = bb[c(2, 4)])
# Plot the choropleth layer
choroLayer(epci_31, var = "nb_rest_10000inhab",
           breaks = bks, col = cols, border = "grey80", lwd = 0.5,
           legend.pos = "topleft",add = TRUE,
           legend.title.txt = "Number of restaurants\nfor 10,000 inhabitants")
# Plot proportionnal symbols
propSymbolsLayer(epci_31, var="nb_of_rest", col="#440170",border=NA,
                 legend.pos="left", inches=0.4, add = TRUE,
                 legend.title.txt = "Number of restaurants")
# Add a layout layer
layoutLayer(title = "Restaurants", sources = "Insee, 2018",
            author = "Kim & Tim, 2018",
            theme = "green.pal", col = "darkred",
            coltitle = "white", postitle = "center",
            frame = TRUE, scale = 10)
# Add a north (south) arrow
north(pos = "topright", south = TRUE)

With ggplot2:

Solution

library(ggplot2)

map_ggplot <- ggplot() +
  geom_sf(data = fra, colour = "ivory3",
          fill = "ivory") +
  geom_sf(data = epci_31, aes(fill = typo), colour = "grey80") +
  scale_fill_manual(name = "Number of restaurants\nfor 10,000 inhabitants",
                    values = cols, na.value = "#303030")+
   geom_sf(data = epci_31 %>%  st_centroid(),
           aes(size= nb_of_rest), color = "#440154CC", show.legend = 'point')+
  scale_size(name = "Number of restaurants",
             breaks = c(1, 500, 3200),
             range = c(0,18))+
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(epci_31)[c(1,3)],
           ylim = st_bbox(epci_31)[c(2,4)]
  ) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "azure",color=NA)) +
  labs(
      title = "Restaurants",
      caption = "Insee, 2018\nKim & Tim, 2018"
  )

plot(map_ggplot)

With tmap:

TODO Timothée

cartography

ggplot2

tmap

TODO Timothée

Exercise 3 : Interactive maps

1

Load the dataset ‘sir_31’ used previously and map the more than 4,000 restaurants of department 31 with the mapview package. Try using different parameters to customize your map.

Information

For example, you can use the map.types, col.regions, label, color, legend, layer.name, homebutton, lwd … parameters of the mapview function.

Solution

library(mapview)
library(sf)

sir_31 <- readRDS("data/sir_31.rds")
mapview(sir_31, map.types = "OpenStreetMap",
        col.regions = "#940000",
        label = paste(sir_31$L2_NORMALISEE, sir_31$NOMEN_LONG, sep = " - "),
        color = "white", legend = TRUE, layer.name = "Restaurants in SIRENE",
        homebutton = FALSE, lwd = 0.5)

** TODO To continue Comeetie **

reproducibility

sessionInfo()

R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.7.0
LAPACK: /usr/lib/lapack/liblapack.so.3.7.0

locale:
 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C               LC_TIME=fr_FR.UTF-8       
 [4] LC_COLLATE=fr_FR.UTF-8     LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] mapview_2.6.3     ggplot2_3.1.0     bindrcpp_0.2.2    cartography_2.2.1
[5] dplyr_0.7.8       sf_0.7-2          knitr_1.21        stringr_1.3.1    

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5   xfun_0.4           purrr_0.3.0        lattice_0.20-38   
 [5] colorspace_1.4-0   viridisLite_0.3.0  stats4_3.5.2       htmltools_0.3.6   
 [9] yaml_2.2.0         base64enc_0.1-3    rlang_0.3.1        e1071_1.7-0.1     
[13] pillar_1.3.1       later_0.7.5        glue_1.3.0         withr_2.1.2       
[17] DBI_1.0.0          RColorBrewer_1.1-2 sp_1.3-1           bindr_0.1.1       
[21] plyr_1.8.4         rgeos_0.4-2        munsell_0.5.0      gtable_0.2.0      
[25] unilur_0.4.0.9000  raster_2.8-19      htmlwidgets_1.3    codetools_0.2-15  
[29] evaluate_0.12      httpuv_1.4.5.1     crosstalk_1.0.0    class_7.3-14      
[33] Rcpp_1.0.0         xtable_1.8-3       satellite_1.0.1    scales_1.0.0      
[37] promises_1.0.1     classInt_0.3-1     jsonlite_1.6       webshot_0.5.1     
[41] leaflet_2.0.2      mime_0.6           png_0.1-7          digest_0.6.18     
[45] stringi_1.2.4      shiny_1.2.0        grid_3.5.2         tools_3.5.2       
[49] magrittr_1.5       lazyeval_0.2.1     tibble_2.0.1       crayon_1.3.4      
[53] pkgconfig_2.0.2    assertthat_0.2.0   rmarkdown_1.11     R6_2.3.0          
[57] units_0.6-2        compiler_3.5.2

In French, IRIS is an acronym of ‘aggregated units for statistical information’. Their target sizes are 2000 residents per basic unit.↩

Let’s Practice Spatial data and cartography with R

Kim A. Etienne C. Timothée G.

SatRday Paris February 2019

Exercise 1 : Manipulate sf objects and associated data.frames

1

Clues

Solution

2

Solution

Solution

3

Solution

Solution

4

Solution

Solution

5

Clues

Solution

Solution

6

Information

Clues

Solution

Exercise 2 : Static maps

1

Information

Solution

2

Solution

Solution

cartography

ggplot2

tmap

Exercise 3 : Interactive maps

1

Information

Solution